#include "cubaturs.h"

#define EM_E0 GSL_CONST_MKSA_VACUUM_PERMITTIVITY //8.854187817e-12
#define EM_M0 GSL_CONST_MKSA_VACUUM_PERMEABILITY //1.25663706144e-6	//electric&magnetic constants
#define EM_LS GSL_CONST_MKSA_SPEED_OF_LIGHT //2.99792458e8		//скорость света

int j;

Cubaturs::Cubaturs()
{
    N=100;   //number of steps for cubaturs
    mI=1; nI=1;
    waveType=2; //2: H-type; if 1: E
    e1=1; e2=1; m1=1; m2=1;
}

Cubaturs::Cubaturs(QString fileName, int method)
{
    QFile readFile(fileName);
    if (!readFile.open(QIODevice::ReadOnly))
        return;
    QTextStream stream(&readFile);
    N=100;   //number of steps for cubaturs
    //a=stream.readLine().toDouble();
    a=stream.readLine().split("=").takeLast().toDouble();
    b=stream.readLine().split("=").takeLast().toDouble();
    if(method==2)
    {
        R=stream.readLine().split("=").takeLast().toDouble();
        fi=stream.readLine().split("=").takeLast().toDouble();
        Z=stream.readLine().split("=").takeLast().toDouble();
        Dist=Z;
    }
    else
    {
        R=stream.readLine().split("=").takeLast().toDouble();
        teta=stream.readLine().split("=").takeLast().toDouble();
        fi=stream.readLine().split("=").takeLast().toDouble();
        Dist=Z;
    }
    e1=stream.readLine().split("=").takeLast().toDouble();
    e2=stream.readLine().split("=").takeLast().toDouble();
    //fi*=M_PI/180;
    //teta*=M_PI/180;
    double length=0.0032;
    k=2*M_PI/length;
    w=k*EM_LS;
    f0=w/2/M_PI;
    m1=1; m2=1;
    alfa();
    beta=sqrt(pow(k,2)*e1*m1-pow(alfa(),2));
    alfa2=1I*sqrt(pow(beta,2)-pow(k,2)*e2*m2);
    I1b=gsl_sf_bessel_I1(__imag__(alfa2)*b);
    I0b=gsl_sf_bessel_I0(__imag__(alfa2)*b);
    K1b=gsl_sf_bessel_K1(__imag__(alfa2)*b);
    K0b=gsl_sf_bessel_K0(__imag__(alfa2)*b);
    G1=G(1);
    G2=G(2);
}

void Cubaturs::calc(int area, int method)
{
    double ro0=0, fis0=0, ro, fis, rr;
    //r0, fis - Переменные интегрирования
    //ro0, fis0 - узлы сетки
    //rr - r при данном ro
    double dq=a/N, dp=2*M_PI/N;	//Приращение по fi и по ro соответственно
    double bb[2]={ 0.2113248654, 0.7886751346 };
    //bb - контейнер кубатуры, bb=[1+beta,1-beta]
    __complex__ double ii1[9], ia=0I, ib=0I, ic=0I, ri,ss[9];
    //ss[0-2] - iiErR,iiErT,iiErF,
    //ss[3-5] - iiEgR,iiEgT,iiEgF,
    //ss[6-8] - iiEjR,iiEjT,iiEjF;
    //ii1 - внутренняя интегральная сумма
    //ia,ib,ic - вспомогательные суммы
    //ri - выражение (1/r+ik)

    fi*=M_PI/180;
    teta*=M_PI/180;
    for(int j=0;j<9; j++) {	//Инициализация
        RES[j]=0I; ss[j]=0I; ii1[j]=0; }
    ee=e1*EM_E0; mm=m1*EM_M0;
    if(area==2)
    {   // for 2-layer second lay
        ro0=a;
        dq=(b-a)/N;
        ee=e2*EM_E0; mm=m2*EM_M0;
    }
    for(int p=0;p<N;p++)
    {
        for(int q=0;q<N;q++)
        {
            for(int s1=0;s1<2;s1++)
            {
                fis=fis0+bb[s1]*dp;
                for(int s2=0;s2<2;s2++)
                {
                    ro=ro0+bb[s2]*dq;
                    getTangPole(area,ro,fis);
                    switch (method)
                    {
                    case 0 :
                        //Численное интегрирование по строгому алгоритму
                        rr=sqrt(pow(R,2)+pow(ro,2)-2*R*ro*cosg(fis));
                        ri=(1/rr+1I*k);
                        //ErR
                        ii1[0]=Ero*ro/rr*ri;
                        ii1[0]*=ro*cexp(rr);
                        ss[0]+=ii1[0];
                        //ErT
                        ii1[1]=1/rr*ri*(R*(Ero*cos(fi-fis)+Efi*sin(fi-fis))-ro*Ero*sin(teta));
                        ii1[1]*=ro*cexp(rr);
                        ss[1]+=ii1[1];
                        //ErF
                        ii1[2]=ri/rr*(Ero*cosb(fis)-Efi*cosd(fis));
                        ii1[2]*=ro*cexp(rr);
                        ss[2]+=ii1[2];
                        //Public multiplier for grad_div
                        ic=3/rr*ri-pow(k,2);
                        ic*=(R*(Hro*cosa(fis)-Hfi*cosg(fis))+ro*Hfi)/pow(rr,2);
                        //EgR
                        ia=(R-ro*cosg(fis))*ic;
                        ib=ri/rr*(Hro*cosa(fis)-Hfi*cosg(fis));
                        ii1[3]=ia-ib;
                        ii1[3]*=ro*cexp(rr);
                        ss[3]+=ii1[3];
                        //EgT
                        ia=ro*cosd(fis)*ic;
                        ib=ri/rr*(Hro*cosb(fis)-Hfi*cosd(fis));
                        ii1[4]=ia+ib;
                        ii1[4]*=ro*cexp(rr);
                        ss[4]+=ii1[4];
                        //EgF
                        ia=ro*sin(fi-fis)*ic;
                        ib=ri/rr*(Hro*cos(fi-fis)+Hfi*sin(fi-fis));
                        ii1[5]=ia-ib;
                        ii1[5]*=ro*cexp(rr);
                        ss[5]+=ii1[5];
                        //EjR
                        ii1[6]=Hro*cosa(fis)-Hfi*cosg(fis);
                        ii1[6]*=ro*cexp(rr);
                        ss[6]+=ii1[6];
                        //EjT
                        ii1[7]=Hro*cosb(fis)-Hfi*cosd(fis);
                        ii1[7]*=ro*cexp(rr);
                        ss[7]+=ii1[7];
                        //EjF
                        ii1[8]=Hro*cos(fi-fis)+Hfi*sin(fi-fis);
                        ii1[8]*=ro*cexp(rr);
                        ss[8]+=ii1[8];
                        break;
                    case 1 :
                        rr=sqrt(pow(R,2)+pow(ro,2)-2*R*ro*cosg(fis));
                        ri=(1/rr+1I*k);
                        //Строгий алгоритм для дальней зоны
                        // What is it? di=(1/(R-ro*cosg(fis))+1I*k);
                        //ErR
                        ii1[0]=Ero*ro;
                        ii1[0]*=ro*cexp0(ro,fis);
                        ss[0]+=ii1[0];
                        //ErT
                        ii1[1]=ro*(Ero*cosa(fis)+Efi*cosg(fis))*sin(fi-fis);
                        ii1[1]*=ro*cexp0(ro,fis);
                        ss[1]+=ii1[1];
                        //ErF
                        ii1[2]=ro*(Ero*cosa(fis)+Efi*cosg(fis))*cosd(fis);
                        ii1[2]*=ro*cexp0(ro,fis);
                        ss[2]+=ii1[2];
                        //EgR
                        ii1[3]=(Hro*cosa(fis)-Hfi*cosg(fis))*(2-1I*k*ro*cosg(fis)-1I*k*ro*Hfi);
                        ss[3]+=ii1[3];
                        //EgT
                        ia=(Hro*cosb(fis)-Hfi*cosd(fis))*(2+pow(k*ro,2)*cosa(fis)*cosg(fis));
                        ib=pow(k*ro,2)*Hfi*cosd(fis);
                        ii1[4]=ia-ib;
                        ii1[4]*=ro*cexp0(ro,fis);
                        ss[4]+=ii1[4];
                        //EgF
                        ia=1I*k*ro*cosa(fis)*sin(teta)*Hro+pow(k*ro,2)*cosb(fis)*cos(teta)*Hfi;
                        ib=(Hro*cos(fi-fis)+Hfi*sin(fi-fis))*(2-pow(k*ro,2)*cosa(fis)*sin(teta));
                        ii1[5]=ib-ia;
                        ii1[5]*=ro*cexp0(ro,fis);
                        ss[5]+=ii1[5];
                        //EjR
                        ii1[6]=Hro*cosa(fis)-Hfi*cosg(fis);
                        ii1[6]*=ro*cexp0(ro,fis);
                        ss[6]+=ii1[6];
                        //EjT
                        ii1[7]=Hro*cosb(fis)-Hfi*cosd(fis);
                        ii1[7]*=ro*cexp0(ro,fis);
                        ss[7]+=ii1[7];
                        //EjF
                        ii1[8]=Hro*cos(fi-fis)+Hfi*sin(fi-fis);
                        ii1[8]*=ro*cexp0(ro,fis);
                        ss[8]+=ii1[8];
                        break;
                    case 2 :
                        rr=sqrt(pow(R,2)+pow(ro,2)+pow(Z,2)-2*R*ro*cos(fi-fis));
                        ri=(1/rr+1I*k);
                        //Строгий алгоритм для цилиндрической системы координат
                        //ErR
                        ii1[0]=(Ero*cos(fi-fis)+Efi*sin(fi-fis))*ri/rr;
                        ii1[0]*=ro*cexp(rr);
                        ss[0]+=ii1[0];
                        //ErF
                        ii1[1]=(Ero*sin(fi-fis)-Efi*cos(fi-fis))*ri/rr;
                        ii1[1]*=ro*cexp(rr);
                        ss[1]+=ii1[1];
                        //ErZ
                        ii1[2]=(ro*Ero-R*(Ero*cos(fi-fis)+Efi*sin(fi-fis)))*ri/rr;
                        ii1[2]*=ro*cexp(rr);
                        ss[2]+=ii1[2];
                        //Multiplier for all grad_div
                        ic=1/pow(rr,2)+ri*(2/rr+1I*k);
                        ic*=(R*(Hro*sin(fi-fis)-Hfi*cos(fi-fis))+ro*Hfi)/pow(rr,2);
                        //EgR
                        ii1[3]=(R-ro*cos(fi-fis))*ic-ri/rr*(Hro*sin(fi-fis)-Hfi*cos(fi-fis));
                        ii1[3]*=ro*cexp(rr);
                        ss[3]+=ii1[3];
                        //EgT
                        ii1[4]=ro*sin(fi-fis)*ic-ri/rr*(Hro*cos(fi-fis)+Hfi*sin(fi-fis));
                        ii1[4]*=ro*cexp(rr);
                        ss[4]+=ii1[4];
                        //EgF
                        ii1[5]=ic;
                        ii1[5]*=ro*cexp(rr);
                        ss[5]+=ii1[5];
                        //EjR
                        ii1[6]=Hro*sin(fi-fis)-Hfi*cos(fi-fis);
                        ii1[6]*=ro*cexp(rr);
                        ss[6]+=ii1[6];
                        //EjT
                        ii1[7]=Hro*cos(fi-fis)+Hfi*sin(fi-fis);
                        ii1[7]*=ro*cexp(rr);
                        ss[7]+=ii1[7];
                        //EjF
                        ii1[8]*=1;
                        ii1[8]*=ro*cexp(rr);
                        ss[8]+=ii1[8];
                        break;
                    }
                }
            }
            for(int j=0; j<9; j++) RES[j]+=ss[j];	//Прибавление результата по очередной кубатуре
            for(int j=0;j<9; j++) ss[j]=0I;	//инициализация нового узла
            ro0+=dq;		//инициализация нового узла
        }
        fis0+=dp; //инициализация новой строки кубатур
        ro0=0;	//возвращение каретки
    }
    for(int j=0; j<9; j++) RES[j]*=dp*dq/4;	//ending of Integer
    switch (method)
    {
    case 0 :
        //Численное интегрирование по строгому алгоритму
        RES[0]*=(cos(teta)/4/M_PI);
        RES[1]/=(4*M_PI);
        RES[2]*=(-1*R/4/M_PI);
        RES[3]/=(4*M_PI*1I*w*ee);
        RES[4]/=(-4*M_PI*1I*w*ee);
        RES[5]/=(4*M_PI*1I*w*ee);
        RES[6]*=(-1)*1I*w*mm/4/M_PI;
        RES[7]*=(-1)*1I*w*mm/4/M_PI;
        RES[8]*=(-1)*1I*w*mm/4/M_PI;
        break;
    case 1 :
        //Строгий алгоритм для дальней зоны
        RES[0]*=cexp(R)*1I*k*cos(teta)/(4*M_PI*R);
        RES[1]*=cexp(R)*1I*k/(-4*M_PI*R);
        RES[2]*=cexp(R)*1I*k/(-4*M_PI*R);
        RES[3]/=(4*M_PI*R*1I*w*ee);
        RES[4]*=cexp(R)/(4*M_PI*R*R*1I*w*ee);
        RES[5]*=cexp(R)/(4*M_PI*R*R*1I*w*ee);
        RES[6]*=cexp(R)*(-1)*1I*w*mm/4/M_PI;
        RES[7]*=cexp(R)*(-1)*1I*w*mm/4/M_PI;
        RES[8]*=cexp(R)*(-1)*1I*w*mm/4/M_PI;
        break;
    case 2:
        //Строгий алгоритм для цилиндрической системы координат
        RES[0]*=Z/4/M_PI;
        RES[1]*=Z/(-4*M_PI);
        RES[2]/=(4*M_PI);
        RES[3]/=(4*M_PI*1I*w*ee);
        RES[4]/=(4*M_PI*1I*w*ee);
        RES[5]*=Z/(4*M_PI*1I*w*ee);
        RES[6]*=(-1)*1I*w*mm/4/M_PI;
        RES[7]*=(-1)*1I*w*mm/4/M_PI;
        RES[8]*=(-1)*1I*w*mm/4/M_PI;
        break;
    }
}

void Cubaturs::getTangPole(int area, double ro, double fis)
{
    switch (area)
    {
    case 0 :
        Ero=cos(fis)+0I;
        Efi=-1*sin(fis)+0I;
        Hro=sqrt(ee/mm)*sin(fis)+0I;
        Hfi=sqrt(ee/mm)*cos(fis)+0I;
        /*Ero=1+0I;
        Efi=0I;
        Hro=0I;
        Hfi=sqrt(mm/ee)+0I;*/
        break;
    case 1 :
        Ero=Ero1(ro,fis);
        Efi=Efi1(ro,fis);
        Hro=Hro1(ro,fis);
        Hfi=Hfi1(ro,fis);
        break;
    case 2 :
        Ero=Ero2(ro,fis);
        Efi=Efi2(ro,fis);
        Hro=Hro2(ro,fis);
        Hfi=Hfi2(ro,fis);
        break;
    case 3 : //- reserved for main structure
        break;
    case 4 :	// Поле H для равноамплитудного возбуждения
        Hro=sin(fis)+0I;
        Hfi=cos(fis)+0I;
        Ero=sqrt(mm/ee)*cos(fis)+0I;
        Efi=-1*sqrt(mm/ee)*sin(fis)+0I;
        break;
    case 5 :	//радиальноуменьшающееся возбуждение
        Ero=cos(fis)*(1-pow(ro/a,2))+0I;
        Efi=-1*sin(fis)*(1-pow(ro/a,2))+0I;
        Hro=sqrt(ee/mm)*sin(fis)*(1-pow(ro/a,2))+0I;
        Hfi=sqrt(ee/mm)*cos(fis)*(1-pow(ro/a,2))+0I;
        break;
    case 6 :
        alfa1=gsl_sf_bessel_zero_Jnu(mI,nI)/a;
        beta=sqrt(pow(k,2)-pow(alfa1,2));
        Ero=-1*1I*beta/alfa1*besselS(mI,alfa1*ro)*(1+Rln[0]);//*(1+0.4-0.573*1I);    //revert 0.7
        Efi=0I;
        Hro=0I;
        Hfi=-1*1I*w*ee/alfa1*besselS(mI,alfa1*ro)*(1-Rln[0]);//*(1-0.4+0.573*1I);    //revert
        //Ero*=(1-Rln[nI]);
        //Hfi*=(1-Rln[nI]);
        break;
    case 7 :
        alfa1=zerosJs(mI,nI)/a;
        beta=sqrt(pow(k,2)-pow(alfa1,2));
        Ero=0I;
        Efi=1I*w*mm/alfa1*besselS(mI,alfa1*ro)*(1-Rln[0]);//*(1+0.265+1I*0.244);       //revert
        Hro=-1*1I*beta/alfa1*besselS(mI,alfa1*ro)*(1+Rln[0]);//*(1+0.265+1I*0.244);    //0.36 revert
        Hfi=0I;
        break;
    case 8 :
        alfa1=zerosJs(mI,nI)/a;
        Ero=1I*w*mm/pow(alfa1,2)/ro*gsl_sf_bessel_J1(alfa1*ro)*sin(fis)*(1-Rln[0]);//*(1-0.28); //2: 0.28_180
        Efi=1I*w*mm/alfa1*besselS(mI,alfa1*ro)*cos(fis)*(1-Rln[0]);//*(1-0.28);
        Hro=-1*1I*beta/alfa1*besselS(mI,alfa1*ro)*cos(fis)*(1+Rln[0]);//*(1+0.28);
        Hfi=1I*beta/pow(alfa1,2)/ro*gsl_sf_bessel_J1(alfa1*ro)*sin(fis)*(1+Rln[0]);//*(1+0.28);
        break;
    }
}

__complex__ double Cubaturs::Echeck()
{
    __complex__ double Integral;
    double dp, dq;
    double ro0,fis0,ro,fis;
    double bb[2];
    __complex__ double ii1[9],ii2;
    Integral=0I;
    ro0=0; fis0=0;
    dq=a/N; dp=2*M_PI/N;
    bb[0]=0.2113248654;
    bb[1]=0.7886751346;
    ii1[0]=0I; ii2=0I;
    for(int p=0;p<N;p++)
    {
        for(int q=0;q<N;q++)
        {
            for(int s1=0;s1<2;s1++)
            {
                fis=fis0+bb[s1]*dp;
                for(int s2=0;s2<2;s2++)
                {
                    ro=ro0+bb[s2]*dq;
//                    rr=r(ro,fis);
//                    ri=(1/rr+i*k);
//                    r3=ro*R/rr;
                    ii1[0]=1I*ro;
                    ii2+=ii1[0]; //Накопление суммы по кубатуре
                    j++;
                }
            }
            Integral+=ii2;		//Прибавление результата по очередной кубатуре
            ii2=0I; j=0;	//инициализация нового узла
            ro0+=dq;		//инициализация нового узла
        }
        fis0+=dp; //инициализация новой строки кубатур
        ro0=0;	//возвращение каретки
    }
    Integral=Integral*dp*dq/4;
    //EjR[0]*=dp*dq/4*(-1)*i*w*EM_M0*mm/4/M_PI;
    //EjT[0]*=dp*dq/4*(-1)*i*w*EM_M0*mm/4/M_PI;
    //EjF[0]*=dp*dq/4*(-1)*i*w*EM_M0*mm/4/M_PI;
    return Integral;
}
